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We  present  a  method  based  on  the  two-process  model  of  sleep  regula¬ 
tion  for  developing  individualized  biomathematical  models  that  predict 
performance  impairment  for  individuals  subjected  to  total  sleep  loss. 
This  new  method  advances  our  previous  work  in  two  important  ways. 
First,  it  enables  model  customization  to  start  as  soon  as  the  first  per¬ 
formance  measurement  from  an  individual  becomes  available.  This 
was  achieved  by  optimally  combining  the  performance  information  ob¬ 
tained  from  the  individual’s  performance  measurements  with  a  priori 
performance  information  using  a  Bayesian  framework,  while  retaining 
the  strategy  of  transforming  the  nonlinear  optimization  problem  of  find¬ 
ing  the  optimal  estimates  of  the  two-process  model  parameters  into  a 
series  of  linear  optimization  problems.  Second,  by  taking  advantage  of 
the  linear  representation  of  the  two-process  model,  this  new  method 
enables  the  analytical  computation  of  statistically  based  measures  of 
reliability  for  the  model  predictions  in  the  form  of  prediction  intervals. 
Two  distinct  data  sets  were  used  to  evaluate  the  proposed  method. 


Results  using  simulated  data  with  superimposed  white  Gaussian  noise 
showed  that  the  new  method  yielded  50%  to  90%  improvement  in  pa¬ 
rameter-estimate  accuracy  over  the  previous  method.  Moreover,  the 
accuracy  of  the  analytically  computed  prediction  intervals  was  validat¬ 
ed  through  Monte  Carlo  simulations.  Results  for  subjects  representing 
three  sleep-loss  phenotypes  who  participated  in  a  laboratory  study  (82 
h  of  total  sleep  loss)  indicated  that  the  proposed  method  yielded  indi¬ 
vidualized  predictions  that  were  up  to  43%  more  accurate  than  group- 
average  prediction  models  and,  on  average,  10%  more  accurate  than 
individualized  predictions  based  on  our  previous  method. 
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HISTORICALLY,  BIOMATHEMATICAL  PERFORMANCE 
PREDICTION  MODELS  HAVE  FOCUSED  ON  PREDICT¬ 
ING  GROUP- AVERAGE  PERFORMANCE  IMPAIRMENT.1'3 
This  modeling  strategy  necessarily  de-emphasizes  individual 
differences  in  resilience  and  vulnerability  to  sleep  loss  and  con¬ 
tradicts  recent  findings  that  indicate  significant  and  systematic 
differences  among  individuals.4'7  For  example,  on  a  psychomo¬ 
tor  vigilance  test  (PVT)  scale,8  restricted-sleep  laboratory  stud¬ 
ies  show  that  interindividual  differences  account  for  nearly  70% 
of  the  total  variance  in  performance  among  a  group  of  individu¬ 
als,  and  that  the  differences  are  trait-like,  stable,  and  innate  to 
an  individual.5  Hence,  even  if  a  group-average  model  were  ca¬ 
pable  of  accurately  predicting  mean-group  performance,  such 
a  model  would  be  of  limited  benefit  without  knowing  how  this 
translates  into  predictions  at  an  individual  level.9 

Recently,  two  methods  have  been  proposed  for  developing  in- 
dividual-specific  performance  prediction  models  for  individuals 
subjected  to  total  sleep  loss  with  uncertain  initial  states  (i.e.,  ini¬ 
tial  homeostatic  pressure  to  sleep  and  circadian  phase).1011  Both 
methods  employ  the  two-process  model  of  sleep  regulation12’13 
as  the  underlying  parametric  model  structure  and,  because  of 
the  absence  of  quantitative  individual  biomarkers  of  perfor- 


Submitted  for  publication  January,  2009 
Submitted  in  final  revised  form  March,  2009 
Accepted  for  publication  May,  2009 

Address  correspondence  to:  Jaques  Reifman,  PhD,  Senior  Research 
Scientist,  Bioinformatics  Cell,  Telemedicine  and  Advanced  Technology 
Research  Center,  U.S.  Army  Medical  Research  and  Materiel  Command, 
MCMR-TT,  Bldg  363  Miller  Drive,  Fort  Detrick,  MD  21702;  Tel:  (301)  619- 
7915;  Fax:  (301)619-1983;  E-mail:  jaques.reifman@us.army.mil 

SLEEP,  Vol.  32,  No.  10,  2009 


mance,14 15  use  previously  collected  PVT  data  from  an  individual 
to  customize  the  model  for  that  individual.  These  methods,  which 
involve  the  minimization  of  the  error  between  the  model  fit  and 
a  set  of  performance  measurements  from  a  specific  individual, 
nonetheless  differ  in  the  techniques  used  to  adapt  the  two-pro¬ 
cess  model  parameters  to  an  individual.  The  method  proposed  by 
Van  Dongen  et  al.10  is  conceptually  based  on  the  earlier  work  of 
Olofsen  et  al.,16  with  the  original  exponential  model  substituted 
by  the  two-process  model  of  sleep  regulation.  Using  a  mixed- 
effects  regression  framework,  they  separate  out  inter-  and  intra¬ 
individual  variability  in  previously  recorded  performance  data 
from  a  group  of  individuals  and,  using  the  maximum  likelihood 
principle,  estimate  the  group-average  model  parameters  and 
their  corresponding  variances.  Thereafter,  as  new  PVT  perfor¬ 
mance  measurements  for  an  individual,  not  necessarily  studied 
previously,  become  available,  a  Bayesian  framework  makes  use 
of  the  learned  inter-individual  variability  of  the  group  and  the 
individual’s  performance  data  to  continuously  adapt  the  model 
parameters  (and  the  model)  to  that  individual.  In  contrast,  the 
method  proposed  by  Rajaraman  et  al.11  adapts  the  two-process 
model  parameters  for  an  individual  based  solely  on  the  individ¬ 
ual’s  performance  data,  without  requiring  inter-  and  intra-indi¬ 
vidual  variability  information  from  a  group  of  individuals.  After 
an  initial  set  of  performance  data  for  an  individual  is  collected, 
as  each  new  performance  datum  becomes  available,  the  param¬ 
eters  of  the  two-process  model  are  continuously  adapted  to  that 
specific  individual  by  solving  a  constrained,  linear  least  squares 
(LS)  problem,  in  which  the  constraint  provides  a  trade-off  be¬ 
tween  the  goodness  of  the  model  fit  and  the  requirement  that  the 
solution  follows  the  two-process  model.  If  a  new  performance 
measurement  is  unavailable  at  the  “expected”  sampling  period, 
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both  methods  bypass  model  adaptation  and  estimate  future  per¬ 
formance  levels  using  the  most  recent  parameter  values. 

Both  of  these  methods  have  advantages  and  limitations. 
The  method  by  Van  Dongen  et  al.10  enables  adaptation  of  the 
two-process  model  parameters  using  only  a  few  performance 
measurements  from  an  individual.  Moreover,  based  on  the 
probability  distributions  of  the  model-parameter  estimates, 
the  adopted  Bayesian  framework  allows  for  the  empirical  es¬ 
timation  of  95%  confidence  intervals  around  the  predictions. 
Although  these  empirically  computed  confidence  intervals 
are  not  validated  through  Monte  Carlo  simulations,  and  their 
soundness  has  been  questioned,17  they  provide  a  first  step  in 
addressing  a  desired  capability  of  performance  models  that, 
until  now,  has  been  lacking.  However,  some  limitations  do 
exist:  Because  of  the  nonlinearity  of  the  two-process  perfor¬ 
mance  model  with  respect  to  its  parameters,13  the  LS  solution 
for  finding  the  best  estimates  of  the  parameters  involves  solv¬ 
ing  a  nonlinear  programming  (NLP)  problem.  In  such  prob¬ 
lems,  unless  the  objective  function  (i.e.,  the  LS  equation)  is 
known  to  be  convex,  the  NLP  solution  may  be  a  local  mini¬ 
mum.  Therefore,  the  proposed  Bayesian  framework1016  cannot 
guarantee  unique  (optimal)  estimates  of  the  model  parameters. 
This  implies  that  the  proposed  method  may  negate  a  key  con¬ 
vergence  property  of  adaptive,  or  learning,  algorithms  in  that 
the  model-parameters’  estimates  must  converge  asymptoti¬ 
cally  to  their  underlying  true  (unknown)  values  as  the  number 
of  measurements  increases  and  as  the  amount  of  uncertainty 
in  the  measurements  decreases.1718  Furthermore,  obtaining  the 
optimal  parameter  estimates  through  their  proposed  Bayesian 
framework  is  computationally  expensive,  involving  a  brute- 
force,  grid-search  procedure  in  which  the  computational  cost 
increases  as  the  desired  accuracy  of  the  resulting  parameter 
estimate  increases. 

The  method  proposed  by  Rajaraman  et  al.11  addresses  some 
of  these  limitations.  It  recasts  the  otherwise  nonlinear  opti¬ 
mization  problem  into  a  set  of  linear  optimization  problems, 
whose  solution,  if  it  exists,  is  guaranteed  to  be  unique.  The 
model-parameter  estimates  are  obtained  by  solving  regular¬ 
ized  linear  LS  problems,11  avoiding  the  computationally  cost¬ 
ly  grid-search  optimization.  More  importantly,  this  method 
ensures  that  the  model  parameters  asymptotically  converge 
to  their  true  values  as  the  number  of  measurements  increases 
and  the  amount  of  noise  in  the  observed  measurements  de¬ 
creases  to  zero.  Despite  these  advantages,  Rajaraman  et  al.11 
did  not  address  two  shortcomings.  First,  their  method  requires 
the  accumulation  of  a  minimum  number  of  past  performance 
measurements  of  an  individual  before  model-parameter  esti¬ 
mates  and  performance  predictions  can  be  made.  In  theory, 
at  least  13  measurements  are  required.  However,  in  practice, 
due  to  noise  in  the  measured  performance  data,  a  larger  num¬ 
ber  of  measurements  are  generally  needed.  This  limitation 
precludes  application  of  the  method  to  partial/  chronic  sleep 
restriction  scenarios,  in  which,  because  of  intermittent  sleep/ 
wake  periods,  accumulation  of  large  consecutive  performance 
measurements  is  not  possible.  Second,  they  failed  to  com¬ 
pute  any  measure  of  reliability  of  the  model  predictions.  Such 
measures  would  provide  statistical  error  bounds  within  which 
model  predictions  may  be  trusted  for  a  predefined  coverage 
probability,  say,  95%. 


In  this  paper,  we  extended  our  previous  work  by  proposing 
solutions  to  these  two  shortcomings.  To  overcome  the  require¬ 
ment  that  a  minimum  number  of  performance  data  points  be  col¬ 
lected  from  an  individual  before  the  model  can  be  customized 
and  used  to  make  predictions,  we  employed  a  Bayesian-like  ap¬ 
proach  akin  to  that  of  Van  Dongen. 10  In  this  approach,  we  com¬ 
bined  a  priori  performance  information  with  the  information 
contained  in  the  individual’s  measured  performance  data.  As 
a  result,  the  proposed  method  leveraged  a  priori  performance 
information  and  started  adapting  the  model  parameters  and 
making  predictions  for  a  specific  individual  as  soon  as  the  first 
performance  measurement  from  that  individual  became  avail¬ 
able.  However,  by  retaining  the  strategy  of  transforming  the 
nonlinear  optimization  problem  into  a  series  of  linear  problems, 
whose  solution  is  the  exact  solution  of  the  original  nonlinear 
problem,  the  proposed  method  guaranteed  unique,  optimal  es¬ 
timates  of  the  two-process  model  parameters,  avoiding  a  brute- 
force,  grid-search  procedure  for  computing  the  estimates.  To 
quantitatively  estimate  the  reliability  of  the  model  predictions, 
we  reformulated  the  two-process  model,  using  its  linear  repre¬ 
sentation,11  into  an  autoregressive  (AR)  model  framework,19'22 
which  directly  provides  analytical  expressions  for  computing 
statistically  based  error  bounds  around  the  model  predictions  in 
the  form  of  prediction  intervals  (Pis). 

We  used  two  distinct  data  sets  to  evaluate  the  proposed 
method.  The  first  consisted  of  simulated  performance  data  on  a 
PVT  scale  generated  from  the  two-process  model  with  known 
parameters  and  superimposed  white  Gaussian  noise.  Simulated 
data  allowed  us  to  verify  the  convergence  of  the  model-param¬ 
eter  estimates  to  their  true  values,  analyze  the  trade-off  between 
a  priori  performance  information  and  measured  performance 
data,  and  qualitatively  and  quantitatively  assess  the  Pis.  We 
found  that  the  proposed  method  yielded  parameter  estimates 
that  asymptotically  converged  to  their  true  values  as  the  num¬ 
ber  of  performance  measurements  for  an  individual  increased 
and  the  amount  of  uncertainty  (noise)  in  the  measurements  de¬ 
creased.  Moreover,  the  new  method  yielded  improvements  in 
parameter-estimate  accuracy  of  up  to  90%  over  our  previous 
method.11  In  addition,  the  soundness  of  the  computed  Pis  was 
validated  through  Monte  Carlo  simulations. 

The  second  data  set  consisted  of  PVT  lapses  obtained  from  a 
laboratory  study  of  individuals  subjected  to  82  h  of  total  sleep 
loss.  The  laboratory  data  allowed  testing  of  the  proposed  ap¬ 
proach  within  the  context  of  the  inter-  and  intra-individual 
variability  encountered  in  actual  performance  data  and  making 
comparisons  between  individualized  predictions  and  group- 
average  model  predictions.  For  individuals  representing  three 
distinct  sleep-loss  phenotypes  (“vulnerable,”  “average,”  and 
“resilient”),  we  found  that  the  new  method  yielded  individual¬ 
ized  predictions  were  up  to  43%  more  accurate  than  the  group- 
average  model  predictions  and  better  captured  the  circadian  and 
homeostatic  variation  of  the  performance  data.  Moreover,  pre¬ 
diction  accuracy  improved,  on  average,  by  10%  compared  with 
the  previous  method. 

METHODS 

In  this  paper,  we  proposed  a  new  method  that  employed 
the  two-process  model  of  sleep  regulation1213  as  the  underly- 
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ing  model  structure  to  enable  individualized  prediction  of  per¬ 
formance  impairment  due  to  total  sleep  loss.  This  method  is 
based  on  our  previous  work,11  in  which,  by  taldng  advantage 
of  the  linear  representation  of  the  two-process  model,  we  esti¬ 
mate  unique  performance  model  parameters  for  an  individual 
after  collecting  a  minimum  number  of  data  points.  Here,  we 
improved  upon  the  previous  work  by  circumventing  the  re¬ 
quirement  that  a  minimum  number  of  previously  collected  data 
points  be  available  before  model-parameter  estimation  and  per¬ 
formance  predictions  can  commence.  With  this  new  method, 
model  individualization  was  started  as  soon  as  the  first  per¬ 
formance  measurement  became  available.  This  was  achieved 
through  Bayesian  inference,  in  which  a  priori  performance  in¬ 
formation  was  combined  with  the  information  available  from 
the  individual’s  measured  performance  data.  The  a  priori  per¬ 
formance  information  was  obtained  from  two-process  model 
predictions  generated  using  a  priori  model-parameter  values. 
As  each  new  performance  measurement  became  available,  it 
was  added  to  the  individual’s  past  performance  data  and  togeth¬ 
er  used  to  update  and  further  individualize  the  model-parameter 
estimates.  Based  on  the  most  recent  parameter  estimates,  we 
made  predictions  in  accordance  with  a  desired  prediction  ho¬ 
rizon.  We  took  further  advantage  of  the  linear  representation 
of  the  two-process  model  to  reformulate  it  into  an  AR  model 
framework,  which  provides  a  direct  analytical  estimation  of  the 
reliability  of  the  model  predictions  in  terms  of  Pis.19 

Two-Process  Model  of  Sleep  Regulation 

The  two-process  model  of  sleep  regulation  consists  of  two 
separate  processes:12  Process  S  (sleep  homeostasis),  which  is 
dependent  on  sleep/wake  history,  increases  exponentially  with 
wake  time  and  decreases  exponentially  with  sleep/recovery  time 
to  a  basal  value,13  whose  rates  of  increase/decrease  are  individ¬ 
ual-specific,  assumed  to  be  constant,  and  have  unknown  values; 
and  Process  C  (circadian),  which  is  independent  of  sleep/wake 
history  and  represents  a  self-sustaining  oscillator  with  a  24-h 
period.23  The  equations  comprising  the  two-process  model  at 
discrete  sampling  time  index  k  can  be  expressed  as13-24 

S(k)  =  1  -  cxp(-T  / x,.)(l  -  S(k  - 1)),  during  wakefulness,  (1) 

S(k)  =  cxp(-T  / Td)(S(k  -1)),  during  sleep,  (2) 

C(k)  =  sin  - 1)7;  +  ^  (3) 

where  S(k)  denotes  the  value  of  the  sleep  homeostat  at  time  k, 
usually  scaled  between  0  and  1  ;24  C(k)  denotes  a  five-harmonic 
sinusoidal  equation  that  approximates  the  circadian  oscillator 
under  entrained  conditions;25  Ts  represents  the  sampling  period; 
x  and  Tj  represent  the  time  constants  of  Process  S  during  wake¬ 
fulness  and  sleep,  respectively;  x  denotes  the  time  period  of  the 
circadian  oscillator  (~24  h);  ap  where  i  =  1,. .  .,5,  represents  the 
amplitude  of  the  five  harmonics  of  the  circadian  process  (flj  = 
0.97,  a2  =  0.22,  a3  =  0.07,  a4  =  0.03,  and  a5  =  0.001);  and  cp  de¬ 
notes  the  initial  circadian  phase.23 

As  suggested  by  Achermann  and  Borbely,26  we  formulated 
the  temporal  pattern  of  cognitive  performance  as  the  additive 


interaction  of  Processes  S  and  C.  Accordingly,  performance 
P(k)  at  time  k  was  expressed  by 

P(k)  =  aS(k)  +  pC{k),  (4) 

where  a  and  ji  are  real-valued  positive  parameters  that  control 
the  trade-off  between  the  two  processes  of  the  model.  For  total 
sleep  deprivation,  Eq.  (4)  can  be  rewritten  as 

5 

P(k )  =  a  -aS(0)/^1  +^a,  sin[/Vu((k-l)7’s  +  ^)],  (5) 

i=i 

where  y  =  exp (-pT),  with  p  =  l/xr,  to  =  2tt/r  and  a.  =  fia..  Equa¬ 
tion  (5)  is  a  function  of  five  unknown  parameters,  a,  y,  ft,  .S'(0), 
and  cp,  of  which  a,  y,  and  /?  have  been  termed  trait  parameters, 
and  5(0)  and  (p  have  been  termed  state  parameters.1011 

Individual-Specific  Model  Development  Based  on  Performance 

Measurements 


In  our  previous  work,11  we  showed  that  performance  P(k)  in 
Eq.  (5)  is  composed  of  a  linear  combination  of  three  basic  com¬ 
ponents,  a  constant  signal  a,  an  exponential  signal  a5(0)/  2,  and 
five  sinusoidal  signals  ai  sin \ico[(k  - 1)7;  +  <j> )],  where  i  =  1,...,5. 
Moreover,  we  note  that,  if  x(k)  denotes  the  value  of  time  series 
x  at  time  k,  a  shift  operator  Z  can  be  defined  such  that  Z"{x(k )} 
=  x(k+n).  Accordingly,  we  establish  three  linear  operators,  Z-y, 

Z-\  and  ^(Z2  _  2cos {icoTs)Z  + 1), 

1=1 

so  that  when  each  operator  is  individually  applied  to  perfor¬ 
mance  P(k)  in  Eq.  (5),  it  eliminates  one  of  the  three  basic  com¬ 
ponents.  Hence,  by  successively  applying  the  three  operators  in 
any  order  to  P(k)  in  Eq.  (5), 


(z  -  ri  n(Z2  -  2  cos(ia>Ts  )Z  + 1)  (z  - 1) 


P(k)  =  0, 


(6) 


we  eliminate  all  three  basic  components  of  P(k)  and  obtain 
an  identically  zero  equality.  The  expression  within  the  brack¬ 
ets  can  be  expanded  into  a  12,h-order  polynomial  in  Z,  which, 
given  the  definition  of  Z,  represents  a  12th-order  linear,  forward- 
difference  equation.  That  is,  for  any  set  of  values  for  the  five 
parameters  [a,  y,  /?,  5(0),  and  <p ],  P(k)  in  Eq.  (5)  is  a  solution  of 
the  12th-order  difference  equation  represented  by  Eq.  (6). 

To  adapt  these  parameters  and  develop  individual-specific 
models,  our  previous  method  only  makes  use  of  current  and 
past  performance  measurements,  such  as  PVT  lapses,  from  the 
individual  being  modeled.  In  this  formulation,  the  relationship 
between  PVT  performance  measurements y(k),  with  k  =  1,  ..., 
N ,  and  P(k)  in  Eq.  (5)  with  unknown  parameter  values  can  be 
written  as 


y{k)  =  P(k)  +  s{k),  (7) 

where  s(k)  denotes  a  normally  distributed  error  (measurement 
noise)  with  zero  mean  and  variance  d1.  We  uniquely  estimate 
P(k),  or  equivalently,  the  five  parameters  of  the  two-process 
model  in  Eq.  (5)  from  y{k),  by  successively  applying  a  pair  of 
the  three  linear  operators  to  y(k)  and  solving  the  three  resulting 
linear  constrained  LS  problems:  the  first  LS  problem  solves  for 
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y,  the  second  for  a  and  ^(0),  and  the  third  for  /?  and  <p.n  For  discussed  above.  For  example,  y  was  estimated  by  solving  the 

example,  we  estimate  a  unique  value  for  y  by  solving  the  fol-  two-step  constrained  linear  LS  problem 

lowing  two-step  constrained  linear  LS  problem 


.  I 


argmmjargmm 

0<v<1 


|P-y|l  +m  P 


(8) 


where  y  denotes  the  N  x  1  vector  of  performance  measure¬ 
ments  y{k ),  with  k=  1 ,  . . . ,  N\  P  represents  the  corresponding 
A7  x  1  vector  of  the  performance  fit  P{k)\  denotes  the  N-\2 
x  N  matrix  that  represents  the  expression  within  the  brackets 
in  Eq.  (6);  and  A2  is  a  positive  real  number.  To  minimize  Eq. 
(8),  we  first  fix  the  value  of  y  (0  <  y  <  1)  and,  through  con¬ 
strained  LS,11  find  P  that  minimizes  the  expression  within 
the  brackets,  for  a  chosen  l2.  In  this  formulation,  l2  provides 
a  trade-off  for  the  solution  of  P  between  fitting  the  perfor¬ 
mance  measurements  y  (the  first  term  inside  the  brackets) 
and  satisfying  the  constraint  imposed  by  L,  (the  second 
term),27  which  requires  that  P  follows  the  two-process  mod¬ 
el  in  Eq.  (5).  By  setting  l2  to  an  arbitrarily  large  value,  we 
forced  the  solution  of  P  to  follow  the  two-process  model, 
while  simultaneously  fitting  the  measurements  y.  To  avoid 
potential  numerical  problems  with  this  procedure,  we  trans¬ 
formed  Eq.  (8)  to  its  standard  form,28  solved  it  with  l2  set  to 
an  arbitrarily  large  value,  and  transformed  its  solution  back 
to  the  original  problem.  Next,  we  change  the  value  of  y  and 
repeat  this  process  over  the  entire  range  of  possible  values,  0 
<  y  <  1,  to  estimate  the  optimal  y,  which  minimizes  Eq.  (8). 
Using  the  optimal  y  and  the  corresponding  solution  of  P  from 
Eq.  (8),  we  estimate  the  other  four  parameters  of  the  two- 
process  model  by  solving  two  additional  constrained  linear 
LS  problems.11  As  new  performance  measurements  from  the 
individual  being  modeled  become  available,  they  are  com¬ 
bined  with  the  previous  measurements  y  and  together  used  to 
adapt  the  model  parameters  for  that  individual  by  repeating 
the  above  procedure. 

A  limitation  of  this  method  is  the  requirement  that  a  mini¬ 
mum  number  of  past  performance  measurements  be  available 
for  the  individual  before  model-parameter  estimation  and  per¬ 
formance  prediction  can  begin.  In  theory,  at  least  13  measure¬ 
ments  are  required;  however,  in  practice,  due  to  noise  in  the 
measurements  y,  a  larger  number  of  measurements  is  generally 
needed. 


Individual-Specific  Model  Development  Based  on  A  priori 

Knowledge  and  Performance  Measurements 

To  overcome  this  limitation  and  permit  model  individualiza¬ 
tion  and  performance  prediction  as  soon  as  the  first  measured 
performance  datum  becomes  available,  we  modified  the  above 
approach  by  considering  a  priori  performance  information 
in  addition  to  the  information  available  from  the  individual’s 
measured  performance  data.  The  a  priori  performance  infor¬ 
mation  was  obtained  from  performance  data  generated  from 
the  two-process  model  in  Eq.  (5),  with  its  model  parameters 
set  to  a  priori  estimated  values.  This  Bayesian  approach  was 
mathematically  implemented  by  adding  a  term  to  Eq.  (8),  which 
accounted  for  M  a  priori  datay(£)  =  1  -M,  2-M,...fi,  and  solving 
the  augmented  expression  through  the  series  of  optimizations 


argmim  argmin 

(Kycl  ^  P 


||P-y||'  +  A2||p-y|| 


+  A2||LrP 


(9) 


where  y  denotes  the  M  x  1  vector  of  prior  performance  data 
y{k)\p  is  a  positive  real  number;  and  P  represents  the  ( N  +  M)  x 
1  vector  of  the  performance  fit  P  ( k ),  with  k  =  1  -M,  2  -M,...,  N, 
whose  first  M  elements  are  represented  by  P  and  the  remaining 
N  elements  are  represented  by  P.  In  other  words,  P  denotes  the 
concatenation  of  vectors  Pand  P. 

The  first  term  in  Eq.  (9)  quantifies  the  fit  of  the  solution  of  P 
to  the  individual’s  measured  performance  data  y,  whereas  the 
second  term  quantifies  the  fit  to  a  priori  generated  performance 
data  y,  with  degree  p1.  Accordingly,  p2  provides  a  mechanism 
to  obtain  a  solution  of  P  that  provides  a  trade-off  between  our 
trust  in  the  measured  performance  data  and  the  a  priori  perfor¬ 
mance  information.  As  p2  decreases,  we  increase  our  trust  in 
the  measured  data  and  emphasize  the  fit  of  P  to  y,  giving  less 
weight  to  the  prior  information.  This  shift  in  emphasis  increases 
the  variance  of  the  solution  of  P  (i.e.,  it  increases  the  trace  of 
the  covariance  of  P),  since  the  y  measurements  are  noisy.29  In 
the  limit  where  p2  =  0.0,  the  a  priori  performance  information 
is  not  considered,  and  the  last  N  elements  of  the  solution  of 
P  in  Eq.  (9)  are  identical  to  P  obtained  from  solving  Eq.  (8). 
Conversely,  as  p2  increases,  we  increase  our  trust  in  the  prior 
information  and  the  model-parameter  estimates  obtained  from 
the  solution  of  P  converge  to  the  a  priori  selected  values. 

The  optimal  value  for  p2  was  computed  by  the  algorithm  de¬ 
scribed  in  Appendix  A  and  is  a  function  of  a  user-provided  es¬ 
timate  of  the  uncertainty  (i.e.,  noise  level)  cr  in  the  measured 
performance  data  y.  The  computed  p2  was  optimal  in  the  sense 
that,  given  a2,  it  minimized  a  cost  function  that  simultaneously 
accounted  for  the  fit  to  the  measured  data  and  the  fit  to  the  prior 
information.  Accordingly,  as  the  provided  uncertainty  in  the  mea¬ 
sured  data  a2  decreased,  the  optimal  p2  decreased,  accentuating 
the  trust  in  the  measured  data,  and  vice  versa.  The  algorithm  for 
the  optimal  selection  of  p2  endowed  another  desired  feature  of 
Bayesian  estimation  algorithms30  in  that  the  trust  in  the  measure¬ 
ments  (i.e.,  observations)  outweighs  the  trust  in  the  a  priori  infor¬ 
mation  as  the  number  of  measurements  provided  to  the  algorithm 
increases.  As  a  result,  the  optimal  p2  decreased  asymptotically  as 
the  number  of  performance  measurements  N  from  an  individual 
increased.  That  is,  as  more  and  more  performance  measurements 
were  attained  for  an  individual,  the  algorithm  placed  more  and 
more  trust  in  the  measurements,  de-emphasizing  the  prior  infor¬ 
mation.  The  rate  of  decrease  of  p2,  however,  became  faster  or 
slower  as  the  user-provided  uncertainty  a2  in  the  measurements  y 
decreased  or  increased,  respectively. 

We  solved  Eq.  (9)  for  the  optimal  estimate  of  y  following 
steps  similar  to  those  used  in  solving  Eq.  (8)  discussed  above. 
We  performed  this  procedure  over  the  entire  range  0  <  y  <  1  to 
estimate  the  optimal  y.  Next,  using  the  optimal  y  and  the  cor¬ 
responding  solution  of  P,  we  uniquely  estimated  the  other  four 
parameters  of  the  two-process  model  by  solving  the  associated 
constrained  LS  problems.  As  new  performance  measurements 
from  the  individual  became  available,  they  were  combined  with 
the  previous  measurements  y  and  used  together  with  the  prior 
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information  y  to  adapt  the  model-parameter  estimates  and  make 
new  predictions. 

Prediction  Uncertainty  Quantification 

Under  certain  conditions,  individualized  predictions  may  be 
of  limited  value  unless  they  are  accompanied  by  measures  of 
reliability  of  prediction  in  the  form  of  statistically  based  error 
bounds,  such  as  Pis.31  One  way  to  estimate  these  Pis  within 
the  context  of  the  two-process  model  is  to  first  compute  con¬ 
fidence  intervals  for  the  model-parameter  estimates  and  then 
translate  them  into  Pis  through  Eq.  (5).  However,  the  nonlin¬ 
ear  relationship  between  the  two-process  model  and  its  param¬ 
eters13  precludes  the  derivation  of  analytical  functions  that  map 
confidence  intervals  of  the  parameters’  estimates  onto  Pis  for 
the  model  predictions.32  Rather,  Pis  are  often  estimated  through 
computationally  expensive  Monte  Carlo  simulations.33 

Here,  we  analytically  computed  Pis  about  the  model  predic¬ 
tions  by  taking  advantage  of  the  linear  representation  of  the 
two-process  model  in  Eq.  (5).  Using  the  property  that  P(k)  in 
Eq.  (5)  is  a  solution  of  the  12th-order  difference  equation  in 
Eq.  (6), 11  we  considered  P(k)  as  a  12th-order  autoregressive 
(AR)  process  (see  Appendix  B)  in  which,  for  each  time  k,  P(/<) 
is  expressed  as  a  linear  combination  of  previous  performance 
measurements,19,22’34  such  that 

P(k)  =  fjbiP(Jc-i),  (10) 

(=i 

where  b.,  with  i  =  1,. .  .,12,  denotes  the  AR-model  coefficients. 
By  considering  the  correspondence  between  Eqs.  (10)  and  (6) 
[see  Appendix  B],  we  obtained  b.  as  the  coefficients  of  the  12th- 
order  polynomial  given  by  the  expression  within  the  brackets 
in  Eq.  (6).  Finally,  by  taking  advantage  of  the  availability  of  an 
analytic  expression  to  compute  statistically  based  error  bounds 
for  AR-model  predictions,19  we  obtained  the  following  equation 
for  computing  Pis  with  a  coverage  probability  of  100(1-0)  % 

PI„(k)  =  P(k)  ±  Z0JbXW+a\  (11) 

where  0  represents  the  significance  level;  Zg/2  represents  the 
percentage  point  of  a  standard  normal  distribution  with  a  pro¬ 
portion  (9/2  above  it;  b  denotes  the  1x12  vector  of  coefficients 
bp  with  /  =  1,...,12;  P(k)  denotes  the  performance  prediction  at 
time  k  given  previous  measurements  y{k- 1 ),  y(k- 2), . . .  ;E  denotes 
the  covariance  matrix  of  P(k-l),...yP(k-l2),  obtained  by  solving 
Eq.  (9);  and  S2  denotes  the  user-provided  noise-level  estimate  of 
the  performance  measurements  y.  We  computed  P(k)  through 
Eq.  (10),  where  P(k-i),  i  =  1,...,12,  were  obtained  from  the  solu¬ 
tion  of  P  in  Eq.  (9).  To  compute  PI ]mul;]%(k-n),  for  n  >  1,  we 
recursively  applied  Eq.  (10)  to  generate  predictions  at  k  +  m, 
with  m  =  1,2 followed  by  Eq.  (11)  for  computing  the  cor¬ 
responding  Pis. 

RESULTS 

We  employed  two  data  sets  to  validate  our  proposed  ap¬ 
proach.  The  first  consisted  of  simulated  performance  data  ob¬ 
tained  from  running  the  two-process  performance  model  in  Eq. 


(5)  with  its  parameters  set  to  selected  values  and  superimposing 
random  noise  to  emulate  the  uncertainty  present  in  measured 
performance  data.  Simulated  data  enabled  verification  of  key 
convergence  properties  of  the  proposed  approach,  sensitivity 
analyses  of  the  trade-off  between  prior  information  and  mea¬ 
sured  data,  and  quantitative  and  qualitative  assessment  of  the 
analytically  computed  Pis. 

The  second  data  set  consisted  of  PVT  lapse  measurements 
from  nine  healthy  individuals  who  had  participated  in  a  82-h  to¬ 
tal  sleep  deprivation  study  under  laboratory  conditions.11,35  The 
laboratory  data  set  allowed  testing  of  the  proposed  approach 
within  the  context  of  the  inter-  and  intra-individual  variabil¬ 
ity  that  characterizes  measured  performance  data  and  compar¬ 
ing  the  individualized  predictions  versus  those  obtained  with 
group-average  prediction  models. 

Simulated  Data  Set 

We  generated  simulated  data  sets  by  running  the  two-process 
performance  model  in  Eq.  (5)  with  known  parameters  and  su¬ 
perimposing  white  (i.e.,  uncorrelated)  Gaussian  noise  to  emu¬ 
late  the  variability  observed  in  measured  performance  data.  By 
sampling  from  the  group-average  distribution  of  parameters  es¬ 
timated  by  Van  Dongen  et  al.,10  we  set  the  three  trait  parameters 
in  Eq.  (5)  to  a  =  30.30,  p  =  6.35,  and  p  =  0.03  h  1  (y  =  0.94), 
whereas  by  sampling  from  uniform  distributions  with  support 
[0  1]  and  [0  24],  respectively,  we  set  the  state  parameters  S(0) 
=  0.82  and  <p  =  6  h.  In  each  simulation,  performance  data  were 
generated  for  82  h  of  total  sleep  deprivation  and  sampled  once 
every  2  h  to  emulate  the  sampling  rate  of  performance  mea¬ 
surements  in  typical  sleep  studies.  Data  from  each  simulation 
were  superimposed  with  white  noise  sampled  from  a  Gaussian 
distribution  with  zero  mean  and  one  of  three  levels  of  variance: 
16,  4,  or  1.  These  variances  at  82  h  of  wakefulness  correspond¬ 
ed  to  a  signal-to-noise  ratio  (SNR)  of  3.70  (59.10/16),  14.77 
(59.10/4),  and  59.10  (59.10/1),  respectively,  where  59.10  was 
the  variance  of  the  performance  data  generated  by  running  Eq. 
(5)  with  the  abovementioned  parameters  for  82  h.  The  SNR  is 
defined  as  the  ratio  of  the  variance  of  a  signal  to  the  variance  of 
the  white  noise  corrupting  that  signal.  In  each  simulation,  the  a 
priori  values  of  the  trait  parameters,  a  r  =  29.70,  [ipr  =  4.30,  and 
ppr  =  0.03  h  1  (y  r  =  0.94),  were  chosen  to  match  Van  Dongen’s 
estimates  of  their  group-average  model  parameters,10  whereas 
the  a  priori  values  of  the  state  parameters,  S(0)pr  =  0.92  and 
<ppr  =  12.6  h,  were  randomly  chosen  by  sampling  from  uniform 
distributions,  as  above.  These  same  parameters  were  used  in 
the  group-average  model  predictions  for  all  simulations.  Based 
on  these  values,  we  employed  Eq.  (5)  to  generate  13  prior  per¬ 
formance  data  points.  Thirteen  is  the  minimum  number  of  data 
points  required  to  recover  the  two-process  model-parameter 
values,  since  the  model  is  the  solution  of  the  12th-order  dif¬ 
ference  equation  [i.e.,  Eq.  (6)].  Although  a  greater  amount  of 
previously  collected  data  may  be  used,  we  empirically  found 
that  13  data  points  provided  the  fastest  rate  of  convergence  of 
the  model-parameter  estimates  to  their  true  values. 

To  make  the  first  prediction,  the  13  a  priori  data  were  com¬ 
bined  with  the  first  performance  measurement  j(l)  to  estimate 
the  homeostatic-rate  parameter  p,  or  equivalently  y,  by  mini¬ 
mizing  Eq.  (9).  We  set  l2  to  21024  (the  largest-possible  double 
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Figure  1 — Simulated  performance  data,  normalized  to  a  psychomotor  vigilance  task  scale,  for  82  h  of  total  sleep  deprivation  with  a  signal-to- 
noise  ratio  of  14.77.  The  solid  circles  show  the  simulated  performance  lapses  every  2  h.  The  dotted,  dash-dotted,  and  dashed  lines  represent  2-,  6-, 
and  10-h-ahead  predictions,  respectively,  based  on  individualized  models  with  prior  parameter  values  as  those  used  in  the  group-average  model 
and  noise-level  estimate  set  to  the  correct  value  (a2  =  4).  The  solid  line  represents  the  group-average  model  predictions.  We  used  the  root  mean 
squared  error  between  the  predicted  and  the  simulated  performance  data,  calculated  from  10  through  82  h,  as  a  comparative  metric. 


precision  floating  point  number)  and,  given  the  user-provided 
uncertainty  <r  in  j(l),  estimated  the  optimal  p2  as  discussed 
in  the  Methods  Section.  The  remaining  four  parameters  of  the 
two-process  model  [a,  5(0),  fi,  and  <p  ],  were  estimated  by  solv¬ 
ing  the  associated  LS  problems,11  with  P  taken  as  the  solution  of 
Eq.  (9).  Thereafter,  as  new  performance  measurements  became 
available,  they  were  added  to  previous  measurements,  one  at  a 
time,  and,  by  following  the  procedure  above,  the  model  param¬ 
eters  were  adapted  to  take  into  account  the  entire  set  of  mea¬ 
surements  up  to  the  last  datum.  Using  the  most  recent  estimates 
of  the  parameters,  we  made  performance  predictions  for  a  de¬ 
sired  prediction  horizon.  For  example,  a  6-h-ahead  prediction 
at  time  index  k  was  made  using  the  model-parameter  estimates 
obtained  at  time  index  k- 3. 

Figure  1  shows  a  representative  realization  of  the  perfor¬ 
mance  profile  of  a  simulated  individual,  with  SNR  =14.77  at  82 
h  of  wakefulness  (i.e.,  a2  =  4),  and  individualized  performance 
predictions  for  three  selected  prediction  horizons  of  2,  6,  and  10 
h.  The  performance  units  were  normalized  to  emulate  results  of 
typical  10-min  PVT  laboratory  studies.  The  root  mean  squared 
error  (RMSE)  was  used  as  a  metric  to  quantify  the  accuracy  of 
the  predictions  (the  smaller  the  RMSE  between  the  predicted 
and  the  true  underlying  performance  signal,  the  better  is  the 
resulting  prediction).19  For  consistency,  RMSEs  were  computed 
across  overlapping  prediction  time  intervals  spanning  from  10 
to  82  h.  As  shown  in  F igure  1 ,  the  individualized  model  was  con¬ 
sistently  more  accurate  than  the  group-average  model,  yielding 
a  two-  to  three-fold  reduction  in  prediction  errors.  Moreover, 
we  note  that  the  group-average  model  predictions  (solid  line) 
were  out  of  phase  with  the  simulated  performance  profile  (solid 
circles).  As  expected,  the  RMSEs  increased  monotonically  with 
increasing  prediction  horizons.  This  should  be  expected  because 
the  model  parameters  were  not  updated  over  the  prediction  ho¬ 
rizon,  preventing  longer-horizon  predictions  from  employing 
intermediate  measurements  to  adapt  the  model  parameters  and 
improve  the  predictions.  Flowever,  it  should  be  noted  that  uti¬ 


lization  of  intermediate  measurements  with  low  SNR  could  re¬ 
sult  in  incorrect  model-parameter  estimates  that,  while  yielding 
a  good  fit  to  the  measurements,  result  in  poorer  predictions. 18  As 
a  result,  low  SNR  measurements  could  yield  predictions  over  a 
longer  horizon  with  smaller  RMSEs  than  ones  over  shorter  ho¬ 
rizons.  Nonetheless,  as  the  number  of  measurements  increases, 
the  model  parameters  estimates  should  converge  to  their  true 
(unknown)  values,  and  the  RMSEs  should  be  approximately  the 
same  regardless  of  the  prediction  horizon.11 

For  the  simulation  in  Figure  1,  Figure  2  shows  the  param¬ 
eter  estimates  (solid  squares),  their  true  values  (dashed  lines), 
and  a  priori  values  (solid  lines)  of  the  five  model  parameters 
as  a  function  of  the  number  of  available  performance  measure¬ 
ments,  starting  with  the  first  PVT  measurement  at  time  zero.  At 
first,  the  parameter  estimates  oscillate  due  to  the  low  SNR  in  the 
simulated  performance  data  during  this  initial  period.  Later,  as 
additional  measurements  became  available  and  the  variance  in 
the  data  increased,  achieving  a  SNR  =  14.77  at  82  h  of  wakeful¬ 
ness,  each  of  the  parameters  asymptotically  converged  to  their 
true  values,  attaining  95%  accuracy  after  ~50  h  of  wakeful¬ 
ness. 

To  verify  the  ability  of  the  algorithm  to  consistently  improve 
the  estimation  of  the  model  parameters  as  a  function  of  data 
uncertainty  <r  and  the  number  of  performance  measurements 
N,  we  performed  three  sets  of  simulations.  Each  consisted  of 
100  trials,  in  which  each  trial  in  a  set  was  superimposed  with 
white  noise  sampled  from  one  of  the  three  Gaussian  distribu¬ 
tions  with  mean  zero  and  variances  that,  at  82  h  of  wakefulness, 
corresponded  to  SNRs  of  59. 10  (a2  =  1),  14.77  (o2  =  4),  and  3.70 
(cr2  =  16).  In  this  way,  we  observed  the  impact  of  each  noise 
level  on  parameter  estimation.  Table  1  shows  the  mean,  the 
variance,  and  the  mean  squared  error  (MSE),  i.e.,  the  square  of 
the  bias  plus  variance,  for  each  of  the  five  parameter  estimates 
over  the  100  trials  for  the  three  sets  of  simulations.  The  MSE, 
which  accounts  for  both  estimation  accuracy  (bias)  and  estima¬ 
tion  precision  (variance),  is  the  most  statistically  meaningful 
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Figure  2 — Estimates  (solid  squares),  the  true  value  (dashed  lines),  and  the  prior  values  (solid  lines)  of  the  five  parameters  of  the  two-process 
model.  These  correspond  to  the  simulated  performance  data  in  Figure  1  and  show  that  the  model-parameter  estimates  asymptotically  ap¬ 
proached  their  true  values  as  they  were  adjusted  every  2  hours,  as  new  measurements  became  available. 


metric,36  and  can  be  quantified  in  simulated  data  where  the  bias 
in  the  estimate  can  be  determined.  The  parameter  estimates  are 
shown  at  four  points  in  time,  at  22,  42,  62,  and  82  h  of  wakeful¬ 
ness,  for  each  of  the  three  SNRs.  For  each  of  the  parameters, 
the  results  consistently  indicated  that,  for  a  fixed  time  point, 
the  MSE  increased  as  the  SNR  decreased  and,  for  a  fixed  SNR, 
the  MSE  decreased  as  the  number  of  available  performance 
measurements  increased.  The  only  exception  was  for  a  at  22  h 
of  wakefulness,  when  the  MSE  for  SNR  =  59.10  (104.44)  was 
larger  than  the  corresponding  value  for  SNR  =  14.77  (98.49). 


This  anomaly  usually  occurs  during  the  early  phases  of  model- 
parameter  estimations  and  when  the  SNR  is  large.  In  this  case, 
the  algorithm  de-emphasizes  the  prior  performance  information 
and  only  trusts  the  measured  data,  which  do  not  contain  suffi¬ 
cient  information  for  model-parameter  estimation.  As  a  result, 
the  estimates  have  large  uncertainty,  occasionally  causing  the 
MSE  to  be  larger  than  expected. 

We  also  compared  the  results  from  our  previous  work" 
with  the  present  results.  For  the  particular  representative  re¬ 
alization  of  the  simulated  performance  profile  in  Figure  1,  we 
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Table  1 — Parameter  Estimates  of  One  Individual  for  Three  Monte  Carlo  Simulations  of  100  Trials,  Each  with  Three  Different  Signal-to-Noise 

Ratios  (59.10,  14.77,  and  3.70) 

a  (lapses) 

P  (h1) 

P  (lapses) 

S(0) 

?(h) 

True  values 

30.30 

0.03 

6.35 

0.82 

6.00 

Signal-to- 

Time 

noise-ratio 

awake  (h) 

59.10 

22 

29.02  (SD  =  10.14) 

0.04  (SD  =  0.03) 

6.39  (SD  =  0.50) 

0.81  (SD  =  0.06) 

6.12  (SD  =  0.33) 

MSE  =  104.44 

MSE  =  0.00 

MSE  =  0.25 

MSE  =  0.00 

MSE  =  0.12 

42 

30.43  (SD  =  3.88) 

0.03  (SD  =  0.01) 

6.25  (SD  =  0.33) 

0.83  (SD  =  0.02) 

6.02  (SD  =  0.22) 

MSE=  15.10 

MSE  =  0.00 

MSE  =  0.12 

MSE  =  0.00 

MSE  =  0.05 

62 

30.29  (SD  =  1.50) 

0.03  (SD  =  0.00) 

6.29  (SD  =  0.27) 

0.83  (SD  =  0.02) 

6.01  (SD  =  0.18) 

MSE  =  2.25 

MSE  =  0.00 

MSE  =  0.08 

MSE  =  0.00 

MSE  =  0.03 

82 

30.20  (SD  =  0.77) 

0.03  (SD  =  0.00) 

6.29  (SD  =  0.24) 

0.83  (SD  =  0.02) 

5.99  (SD  =  0.14) 

MSE  =  0.60 

MSE  =  0.00 

MSE  =  0.06 

MSE  =  0.00 

MSE  =  0.02 

14.77 

22 

28.48  (SD  =  9.76) 

0.05  (SD  =  0.08) 

6.36  (SD  =  0.90) 

0.82  (SD  =  0.09) 

6.19  (SD  =  0.60) 

MSE  =  98.49 

MSE  =  0.01 

MSE  =  0.81 

MSE  =  0.01 

MSE  =  0.40 

42 

30.28  (SD  =  5.91) 

0.03  (SD  =  0.01) 

6.06  (SD  =  0.67) 

0.83  (SD  =  0.04) 

6.07  (SD  =  0.44) 

MSE  =  34.98 

MSE  =  0.00 

MSE  =  0.53 

MSE  =  0.00 

MSE  =  0.20 

62 

30.26  (SD  =  2.63) 

0.03  (SD  =  0.01) 

6.20  (SD  =  0.54) 

0.83  (SD  =  0.03) 

6.03  (SD  =  0.35) 

MSE  =  6.93 

MSE  =  0.00 

MSE  =  0.32 

MSE  =  0.00 

MSE  =  0.13 

82 

30.14  (SD  =  1.44) 

0.03  (SD  =  0.00) 

6.19  (SD  =  0.47) 

0.83  (SD  =  0.03) 

5.99  (SD  =  0.28) 

MSE  =  2.11 

MSE  =  0.00 

MSE  =  0.25 

MSE  =  0.00 

MSE  =  0.08 

3.70 

22 

26.63  (SD  =  10.19) 

0.06  (SD  =  0.05) 

6.21  (SD=  1.69) 

0.83  (SD  =  0.15) 

6.41  (SD  =  1.15) 

MSE  =  117.33 

MSE  =  0.00 

MSE  =  2.87 

MSE  =  0.02 

MSE  =  1.49 

42 

28.66  (SD  =  5.88) 

0.04  (SD  =  0.02) 

5.56  (SD=  1.30) 

0.85  (SD  =  0.09) 

6.28  (SD  =  0.91) 

MSE  =  37.24 

MSE  =  0.00 

MSE  =  2.31 

MSE  =  0.01 

MSE  =  0.90 

62 

29.81  (SD  =  3.44) 

0.03  (SD  =  0.01) 

5.88  (SD  =  1.09) 

0.85  (SD  =  0.07) 

6.18  (SD  =  0.77) 

MSE  =  12.09 

MSE  =  0.00 

MSE  =  1.41 

MSE  =  0.01 

MSE  =  0.62 

82 

29.93  (SD  =  2.21) 

0.03  (SD  =  0.01) 

5.89  (SD  =  0.95) 

0.84  (SD  =  0.06) 

6.05  (SD  =  0.59) 

MSE  =  5.03 

MSE  =  0.00 

MSE=  1.12 

MSE  =  0.00 

MSE  =  0.35 

The  estimates’ 

mean,  standard  deviation  (SD),  and  mean  squared  error  (MSE),  i.e.,  the  square  of  bias  plus  variance 

are  shown  at  four  time 

points:  22, 42,  62,  and  82  h  of  wakefulness. 

observed  significant  prediction  and  parameter-estimate  im¬ 
provements  with  the  current  method.  In  particular,  the  average 
prediction  RMSE  over  the  three  horizons,  computed  between 
52  and  82  h  of  wakefulness,  was  80%  smaller  with  the  current 
method.  The  average  MSEs  over  the  five  parameters  shown 
in  Table  1  at  62  h  of  wakefulness  were  50%,  90%,  and  93% 
smaller  with  the  new  method  when  the  SNR  was  set  at  59.10, 
14.77,  and  3.70,  respectively,  whereas  at  82  h  of  wakefulness, 
the  corresponding  values  were  45%,  50%,  and  89%  smaller 
with  the  new  method.  These  two  time  points  (62  and  82  h  of 
wakefulness)  were  the  only  common  points  to  both  studies. 
The  simulation  results  indicated  that  the  proposed  method’s 
improvements  in  predictions  and  parameter  estimations  were 
more  significant  when  a2  in  the  performance  data  was  high  and 
the  number  of  available  measurements  for  parameter  estima¬ 
tion  was  small. 

Figure  3  shows  the  performance  profile,  the  10-h-ahead  pre¬ 
dictions,  and  their  corresponding  95%  Pis  for  the  simulated  in¬ 
dividual  in  Figure  lwith  SNRs  set  at  3.70  (top),  14.77  (middle), 
and  59.10  (bottom).  The  10-h-ahead  predictions  at  any  sam¬ 
ple  index  k  were  actually  computed  at  sample  index  k- 5  and 
are  illustrated  in  all  three  panels  by  dashed  lines.  The  dotted 
lines  represent  the  corresponding  upper  and  lower  limits  of  the 


analytically  computed  95%  Pis.  To  corroborate  these  estima¬ 
tions,  we  empirically  computed  the  95%  Pis  (dash-dotted  lines) 
through  100  Monte  Carlo  simulations.  The  results  indicated 
that  while  the  analytically  computed  95%  Pis  were  initially  un¬ 
derestimated,  for  each  of  the  three  SNRs,  they  converged  to  the 
corresponding  empirically  calculated  values  as  the  number  of 
performance  measurements  increased.  This  is  attributed  to  the 
Bayesian  aspect  of  our  method,  where  during  the  early  phases 
of  the  model-parameter  estimates  the  algorithm  places  greater 
trust  on  the  prior  performance  information  than  the  measured 
data,  resulting  in  Pis  that  are  narrower  than  expected.  Flowever, 
as  the  number  of  performance  measurements  for  an  individual 
increased,  the  95%  PI  converged  to  the  corresponding  Monte 
Carlo  simulated  values,  covering  almost  entirely  the  individu¬ 
al’s  performance  measurements.  The  results  also  indicated  that, 
as  expected,  the  width  of  the  analytically  computed  95%  Pis 
decreased  from  the  top  to  the  bottom  panels  as  the  SNR  in  the 
simulated  performance  data  increased. 

To  analyze  the  sensitivity  of  the  a  priori  chosen  values  of  the 
two-process  model  parameters  on  the  performance  prediction, 
we  compared  and  contrasted  predictions  for  three  simulated  in¬ 
dividuals,  each  representing  a  different  sleep-loss  phenotype. 
Figure  4  shows  simulated  performance  profiles  (solid  circles), 
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Figure  3 — Ten-hour-ahead  performance  predictions  (dashed  lines)  and  their  corresponding  analytically  computed  95%  prediction  intervals  (dot¬ 
ted  lines)  for  the  simulated  performance  data  shown  in  Figure  1,  with  signal-to-noise  ratio  (SNR)  set  at  3.70  (top),  14.77  (middle),  and  59.10  (bot¬ 
tom).  The  dash-dotted  lines  in  each  of  the  three  panels  represent  the  95%  prediction  intervals  computed  through  100  Monte  Carlo  simulations. 


10-h-ahead  predictions  (dashed  lines),  and  group-average  mod¬ 
el  predictions  (solid  lines)  for  representative  “average”  (top), 
“vulnerable”  (middle),  and  “resilient”  (bottom)  individuals  with 
SNRs  set  to  14.77.  We  computed  the  group-average  predictions 
using  the  same  parameter  values  as  those  used  for  generating  a 
priori  performance  data.  The  figure  shows  that  the  accuracy  in 
the  10-h-ahead  predictions,  measured  in  terms  of  RMSE,  var¬ 
ied  between  individuals,  in  that  the  larger  the  RMSE  between 
the  group-average  model  predictions  and  the  true  underlying 
performance  of  an  individual,  the  larger  the  error  in  the  indi¬ 
vidualized  performance  predictions  for  that  individual.  In  other 
words,  as  expected,  the  closer  the  prior  parameter  values  were 


to  an  individual’s  true  (unknown)  parameters,  the  more  accu¬ 
rate  the  individualized  predictions  were  for  that  individual. 

Except  in  simulations,  one  does  not  know  the  noise  level  d1 
in  performance  data.  Hence,  to  analyze  the  robustness  of  the 
proposed  individualized  prediction  method  against  errors  in 
noise-level  estimates  a2,  we  compared  the  predictions  for  a  sim¬ 
ulated  individual  using  three  different  a2.  Figure  5  shows  three 
10-h-ahead  predictions  for  the  “average”  individual  in  Figure 
4  with  SNR  equal  to  14.77.  In  each  prediction,  the  noise-level 
estimate  a2  was  set  to  one-hundredth  (dotted),  equal  to  (dashed), 
and  100  times  (dash-dotted)  the  true  noise  level  (o2  =  3.38)  used 
in  simulating  performance  profile. 
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Figure  4 — Individualized  model  predictions  for  three  simulated  individuals,  each  representing  three  different  sleep-loss  phenotypes  [average 
sensitivity  to  sleep  loss  (top),  vulnerable  to  sleep  loss  (middle),  and  resilient  to  sleep  loss  (bottom)].  The  solid  circles  in  each  of  the  panels 
represent  the  simulated  performance  data  with  a  signal-to-noise  ratio  of  14.77.  The  dashed  lines  represent  10-h-ahead  predictions  using  the 
group-average  model  parameter  values  as  the  prior  and  using  the  correct  noise  level.  The  solid  lines  in  each  panel  represent  the  group-average 
model  predictions. 


Inherently,  the  proposed  algorithm  followed  a  key  property 
of  learning  algorithms,18  assigning  greater  trust  to  the  measured 
performance  data  as  the  user-provided  noise-level  estimate  (<r2) 
decreased  to  zero,  and  vice  versa.  Figure  5  shows  that,  while 
over-estimation  of  the  noise  level  results  in  predictions  closer  to 
the  group-average  model  predictions,  underestimation  results  in 
predictions  with  RMSEs  equivalent  to  those  obtained  using  the 


correct  noise-level  estimate.  This  was  because  the  two-process 
model  is  composed  of  sinusoidal  and  exponential  components, 
which  do  not  possess  enough  flexibility  to  yield  performance 
signals  that  exactly  fit  a  set  of  noisy  performance  measurements. 
Hence,  underestimation  of  the  noise  level,  which  increased  the 
trust  on  the  measured  performance  data,  did  not  allow  the  algo¬ 
rithm  to  overfit  the  performance  data.  Although,  for  a  particular 
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Figure  5 — Individualized  model  predictions  for  the  average  sleep-sensitivity  phenotype  in  Figure  4  (top).  The  solid  circles  represent  the 
simulated  performance  data  with  a  signal-to-noise  ratio  of  14.77.  The  dotted,  dashed,  and  dash-dotted  lines  represent  the  10-h-ahead  predic¬ 
tions  based  on  noise-level  estimates  or  set  to  one-hundredth,  equal  to,  and  100  times,  respectively,  the  true  noise  level  (a2  =  3.38).  The  prior 
parameter  values  were  the  same  as  those  used  in  the  group-average  model. 


Time  awake  (h) 

Figure  6 — Mean  and  standard  deviation  of  psychomotor  vigilance  task  (PVT)  lapse  measurements  collected  every  2  hours  for  1 1  individuals 
during  82  h  of  total  sleep  deprivation.  The  solid  line  shows  the  performance  predictions  obtained  with  the  group-average  model. 


realization  of  noise  in  the  performance  data,  predictions  with 
under-estimation  of  the  noise  level  may  yield  smaller  RMSEs 
than  those  with  the  correct  noise-level  estimate,  predictions  av¬ 
eraged  over  a  number  of  noise  realizations  should  yield  smaller 
RMSEs  when  <r  is  closer  to  the  true  noise  level  in  the  measured 
data. 

Laboratory  Data  Set 

The  second  data  set  used  to  test  the  proposed  approach 
was  collected  from  a  controlled  laboratory  study  in  which  48 
healthy  adults  were  kept  awake  for  85  consecutive  hours  to 
determine  the  effects  of  fatigue  countermeasures  on  perfor¬ 
mance  and  alertness  during  prolonged  sleep  deprivation.35  In 


the  present  study,  we  tested  a  subset  of  11  subjects  who  had 
not  been  administered  active  fatigue  countermeasures,  and  for 
whom  PVT  measures  were  collected  every  2  h,  for  a  total  of  42 
measurements.  Figure  6  shows  the  mean  performance  profile 
(solid  circles)  and  the  standard  deviation  of  the  11  subjects  over 
the  82  h  of  PVT  data  collection  and  the  group-average  model 
predictions  (solid  lines).  Overall,  the  trend  suggested  that  both 
PVT  lapses  and  variance  (70%  of  which  could  be  attributable 
to  inter-individual  variability5)  increased  over  time.  Moreover, 
the  figure  shows  that  the  group-average  model  does  not  even 
predict  the  mean-group  performance  well. 

For  comparison  purposes,  we  selected  the  same  three  sub¬ 
jects  as  in  the  previous  study,11  each  representing  one  of  three 
different  sleep-loss  phenotypes:  relatively  vulnerable  to  sleep 
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Figure  7 — Individualized  model  predictions  for  three  subjects  with  different  behavioral  phenotypes  [average  sensitivity  to  sleep  loss  (top), 
vulnerable  to  sleep  loss  (middle),  and  resilient  to  sleep  loss  (bottom)],  with  prior  parameter  values  equal  to  those  used  in  the  group-average 
model  and  noise-level  estimate  a1  set  to  77.60.  The  solid  circles  in  each  of  the  panels  represent  the  measured  psychomotor  vigilance  task 
(PVT)  lapses,  measured  every  2  hours.  The  dashed  lines  represent  the  10-h-ahead  predictions,  whereas  the  dotted  lines  represent  the  corre¬ 
sponding  analytically  computed  95%  prediction  intervals.  The  solid  lines  in  each  panel  represent  the  group-average  model  predictions. 


loss,  relatively  average  sensitivity  to  sleep  loss,  and  relatively 
resilient  to  sleep  loss.  The  three  panels  in  Figure  7  show  the 
observed  PVT  lapses  (solid  circles)  for  the  average  (top),  vul¬ 


nerable  (middle),  and  resilient  (bottom)  subjects  along  with 
the  group-average  model  predictions  (solid  lines),  the  10-h- 
ahead  predictions  (dashed  lines),  and  their  corresponding  95% 
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Pis  (dotted  lines).  In  these  calculations,  we  set  the  noise-level 
estimate  a1  =  77.60,  as  in  Van  Dongen.10  For  consistency,  the 
RMSEs  were  computed  between  10  and  82  h,  and  they  indi¬ 
cated  that  the  individualized  predictions  were  up  to  43%  more 
accurate  than  the  group-average  model  predictions  for  all  three 
subjects,  with  the  measurements  falling  almost  completely 
within  the  corresponding  analytically  computed  95%  Pis.  We 
also  compared  the  10-h-ahead  predictions  with  our  previous 
method11  and  found  that,  on  average,  the  new  method  reduced 
the  prediction  error  by  10%. 

DISCUSSION 

In  this  paper,  we  presented  a  new  method  based  on  the  two- 
process  model  of  sleep  regulation  that  enabled  individualized 
predictions  of  performance  impairment  represented  by  PVT 
lapses  for  subjects  exposed  to  total  sleep  loss.  The  method  ad¬ 
vanced  our  previous  work"  in  two  important  ways.  First,  it  al¬ 
lowed  model-parameter  estimation  and  performance  predictions 
as  soon  as  the  first  performance  measurement  became  available. 
This  was  achieved  by  combining  a  priori  performance  informa¬ 
tion  with  the  individual’s  performance  data  through  Bayesian 
inference.  However,  by  retaining  the  strategy  used  in  our  previ¬ 
ous  work,  where  the  nonlinear  optimization  problem  of  find¬ 
ing  the  best  estimates  of  the  two-process  model  parameters  is 
transformed  into  a  series  of  linear  problems,  the  new  method 
guaranteed  unique  estimates  of  the  five  parameters  of  the  two- 
process  model,  avoiding  brute-force,  grid-search  procedures.10 
Another  important  feature  of  the  new  method  was  that  a  priori 
performance  information  was  optimally  combined  with  the  in¬ 
dividual’s  performance  data  as  a  function  of  a  user-provided 
estimate  of  the  uncertainty  (i.e.,  the  noise  level)  in  the  mea¬ 
surements.  As  a  result,  the  algorithm  increased  its  trust  in  the 
a  priori  performance  information  as  the  estimate  of  the  noise 
level  in  the  measurements  increased,  and  vice  versa.  Moreover, 
as  more  and  more  performance  measurements  were  attained  for 
an  individual,  the  algorithm  placed  more  and  more  trust  in  the 
measurements,  de-emphasizing  a  priori  performance  informa¬ 
tion.  The  rate  of  increase  of  the  trust  placed  in  the  measure¬ 
ments,  however,  became  faster  or  slower  as  the  user-provided 
uncertainty  in  the  measurements  decreased  or  increased,  re¬ 
spectively. 

Second,  the  current  work,  for  the  first  time,  provided  sta¬ 
tistically  based  error  bounds  around  the  model  predictions  in 
the  form  of  prediction  intervals.  This  was  achieved  by  taking 
advantage  of  the  linear  representation  of  the  two-process  model 
proposed  in  our  previous  work,"  which  afforded  the  reformula¬ 
tion  of  the  two-process  model  into  an  equivalent  AR  model  and, 
thus,  provided  analytical  expressions  for  estimating  Pis  about 
the  model  predictions,  bypassing  the  need  to  first  estimate  con¬ 
fidence  intervals  of  the  model  parameter  estimates. 10 

Employing  a  simulated  performance  data  set  generated  from 
the  two-process  model  with  selected  parameters  and  with  su¬ 
perimposed  white  Gaussian  noise,  we  found  that  the  param¬ 
eters’  estimates  asymptotically  converged  to  their  true  values 
as  the  number  of  performance  data  points  increased  and  as  the 
amount  of  noise  in  the  performance  data  decreased,  thereby 
confirming  that  the  proposed  method  satisfied  key  convergence 
properties  of  adaptive  algorithms.18  Moreover,  we  observed  that 


performance  predictions  made  over  shorter  horizons  were  more 
accurate  than  predictions  made  over  longer  horizons.  This  ob¬ 
servation  is  because  of  the  more  frequent  use  of  intermediate 
performance  measurements  by  shorter  horizons  in  updating  the 
model-parameter  estimates.  However,  intermediate  measure¬ 
ments  with  low  SNR  could  result  in  inaccurate  models  that, 
while  yielding  a  good  fit  to  the  measurements,  predict  poor¬ 
ly.18  We  also  analyzed  the  analytically  computed  95%  Pis  and 
found,  as  expected,  that  their  width  was  positively  correlated 
with  the  amount  of  noise  in  the  performance  data.  Moreover, 
as  more  performance  measurements  were  available  to  the  algo¬ 
rithm,  the  Pis  converged  to  those  empirically  obtained  through 
Monte  Carlo  trials,  validating  the  accuracy  of  the  Pis  computed 
by  the  new  method.  Employing  PVT  data  from  a  laboratory 
study,  the  proposed  method  yielded  individualized  predictions 
for  three  sleep-loss  phenotypes  that  were  up  to  43%  more  ac¬ 
curate  than  group-average  model  predictions.  Additionally,  the 
corresponding  95%  Pis  almost  entirely  covered  the  entire  set  of 
measurements. 

When  comparing  the  results  of  the  proposed  method  with 
those  from  our  previous  approach  on  the  same  simulated  data, 
we  observed  average  improvements  in  the  accuracy  of  the 
model-parameter  estimates  of  70%  and  in  the  accuracy  of  the 
performance  predictions  of  80%.  In  particular,  we  found  that 
the  improvements  in  parameter  estimation  accuracy  were  more 
pronounced  during  the  early  phases  of  parameter  estimation 
(typically,  before  obtaining  13  measurements)  and  when  the 
SNR  was  low.  This  was  because  the  algorithm  places  more  trust 
on  the  a  priori  information  during  the  early  stages  of  adapta¬ 
tion  (i.e.,  before  the  performance  impairment  trend  could  be 
completely  learned)  and  when  the  estimated  amount  of  noise 
in  the  data  was  large.  Similarly,  on  the  laboratory  data  set,  the 
new  method  yielded  a  10%  average  reduction  in  the  prediction 
error. 

We  found  the  accuracy  of  the  model-parameter  estimates  and 
of  the  performance  predictions  to  be  a  function  of  the  a  priori 
selected  values  for  the  model  parameters  and  the  noise-level 
estimate  of  the  performance  measurements.  We  empiric  ally  ob¬ 
served  that,  as  expected,  the  closer  the  a  priori  parameter  values 
were  to  the  true  values  underlying  an  individual,  the  more  accu¬ 
rate  were  the  predictions  for  that  individual.  Additionally,  over¬ 
estimation  of  the  noise  level  in  the  measurements  resulted  in 
predictions  biased  toward  the  a  priori  information,  whereas  un¬ 
derestimation  of  the  noise  level  did  not  yield  necessarily  worse 
predictions  than  the  ones  based  on  the  correct  value  of  the  noise 
level.  This  was  because  the  two-process  model  of  sleep  regu¬ 
lation  is  composed  of  sinusoidal  and  exponential  components, 
which  prevent  a  perfect  fit  of  the  model  to  non-smooth  signals, 
such  as  additive  white  noise  in  the  performance  measurements. 
This  precludes  model  overfitting  and  resulting  deterioration  of 
model  predictions. 

Despite  the  advances  made  by  the  new  method  for  individu¬ 
alized  performance  prediction  of  sleep-deprived  individuals,  it 
does  have  limitations.  In  particular,  this  method  suffers  from 
the  same  limitations  as  any  Bayesian  approach,  where  a  “good” 
choice  of  a  priori  parameter  values  and  a  “reasonable”  estimates 
of  the  noise-level  in  performance  data  are  key  for  obtaining  opti¬ 
mal  results. 18  Another  limitation  relates  to  the  underlying  assump¬ 
tion  that  measures  of  performance,  such  as  PVT  lapses,  would  be 
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available  on  a  frequent  basis.  This  may  not  be  possible  in  some 
operational  settings,  since  it  may  not  be  practical  to  discontinue 
work-related  tasks  for  administering  perfonnance  tests.  Further¬ 
more,  we  note  that,  although  PVT  lapse  (selected  as  our  predicted 
variable)  is  recognized  as  a  widely  used  and  sensitive  measure  to 
sleep  loss  and  although  PVT  may  be  the  most  practical  test  for 
some  operational  environments,37  it  may  not  accurately  reflect 
real,  work-related  performance  of  individuals. 

Our  future  modeling  efforts  will  focus  on  predicting  perfor¬ 
mance  under  chronic  sleep  restriction  conditions,  where  wake 
and  sleep  episodes  alternate,  precluding  the  availability  of  a 
large  set  of  contiguous  performance  measurements  from  which 
individualized  models  can  be  obtained.  This  difficulty  can 
be  handled  by  our  proposed  method  in  that  model  adaptation 
would  take  place  incrementally  during  wake  episodes,  and  fu¬ 
ture  performance  levels  would  be  predicted  based  on  the  latest 
parameter  values.  In  addition,  we  will  assess  the  performance 
of  the  proposed  approach  for  the  prediction  of  other  outcome 
measures  of  performance,  such  as  the  Karolinska  sleepiness 
scale,38  the  Stanford  sleepiness  scale,39  the  serial  addition/ sub¬ 
traction  task,40  and  the  digit  symbol  substitution  task.41 

Considerable  modeling  efforts  are  still  required  to  fully  ad¬ 
dress  the  set  of  capability  gaps  identified  at  the  Department  of 
Defense-sponsored  Fatigue  and  Performance  Modeling  Work¬ 
shop.1-3  Flowever,  the  work  described  here  provides  significant 
contributions  toward  closing  the  research  gaps  and  developing 
models  that  accurately  predict  cognitive  performance  impair¬ 
ments  due  to  sleep  deprivation  at  an  individual  level. 
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APPENDIX  A 

Optimal  Selection  of  the  Parameter  p2  in  Eq.  (9)  [main  text] 

In  this  appendix,  we  describe  the  method  employed  to  obtain 
an  optimal  value  for  the  parameter  p2  in  Eq.  (9) 


which  simultaneously  accounts  for  the  fit  to  the  measurements 
y  [the  first  tenn  in  Eq.  (A.  1 )]  and  the  fit  to  the  prior  performance 
data  y(the  second  tenn). 

We  obtained  an  optimal  value  for  ir  by  minimizing  the  ex¬ 
pression  (i.e.,  the  cost  function)  within  the  braces  of  the  follow¬ 
ing  equation 


arg  min  j  |p^  —  y||^  +  tr[cov(P^ )]  j,  (A.2) 

where  Vi  denotes  the  last  N  elements  of  the  solution  of  P  in  Eq. 
(A.l),  for  a  specific  p2  and  for  X2  set  to  an  arbitrarily  selected 
large  number  (21024).  To  avoid  potential  numerical  problems 
with  this  procedure,  we  transformed  Eq.  (A.l)  to  its  standard 
form,28  solved  it  for  a  specific  /r  and  for  X2  set  to  an  arbitrarily 
selected  large  number  (21024),  and  transformed  its  solution  back 
to  the  original  problem  [i.e.,  Eq.  (A.l)].  The  first  term  within 
the  braces  in  Eq.  (A.2)  represents  the  fit  of  the  solution  of  Eq. 
(A.l)  to  the  performance  data  measurements  y,  whereas  the  sec¬ 
ond  term  represents  the  trace  of  the  covariance  of  P  A,  which  is 
equivalent  to  the  fit  of  the  solution  of  P  to  the  a  priori- generated 
performance  data  y,  i.e.,  the  second  term  in  Eq.  (A.l).  Noting 
that  P(/!  and  //j  cov(P)  |  are  functions  of  p2  and  X2,  for  which 
closed-form  expressions  exist,28  and  because  Eq.  (A.2)  has  a 
unique  minimum  (see  below),  the  optimal  p2  can  be  obtained 
through  standard  numerical  unconstrained  optimization  tech¬ 
niques,  such  as  Levenberg-Marquardt  or  Gauss-Newton. 


As  p2  increases,  the  first  term  in  Eq.  (A.2)  increases  monoton- 
ically  while  the  second  term  decreases  monotonically.29  There¬ 
fore,  the  cost  function  in  Eq.  (A.2)  has  at  most  one  minimum,42 
which  corresponds  to  the  optimal  p2.  Because  tr[cov(P/(d]  in 
Eq.  (A.2)  is  proportional  to  the  user-provided  uncertainty  <r  in 
the  measurements  y  (i.e.,  the  noise  level  in  y),  the  optimal  p2 
is  a  function  of  a2.  Thereby,  as  a1  increases,  the  optimal  p2  in¬ 
creases,  accentuating  the  trust  in  the  prior  performance  data  in 
Eq.  (A.l),  and  vice  versa. 

APPENDIX  B 

In  this  appendix,  we  show  the  equivalence  between  the  two- 
process  model  of  performance  given  by  Eq.  (5) 

5 

P(k)  =  a-aS{0)yk'  +^a,.  sin[i©((A:-l)7; +^)],  (B.l) 

1=1 

and  the  12th-order  autoregressive  (AR)  model  given  by  Eq.  (10) 
P{k)  =  fjb,P{k-i).  (B.2) 

;=i 

Note  that  for  any  set  of  values  for  the  five  parameters  [a,  y, 
/?,  5(0),  and  cp\,  P(k)  in  Eq.  (B.l)  is  a  solution  to  the  12th-order 
difference  equation  given  by  the  expression  within  the  brackets 
in  Eq.  (6) 


+  l)j(Z- 

where  Z  denotes  a  shift  operator,  such  that  Z"{P(k)}  =  P(k+n). 
By  expanding  the  linear,  forward-difference  operator  in  Eq. 
(B.3)  into  a  12th-degree  polynomial  in  Z  and  applying  it  to  P(k) 
in  Eq.  (B.l),  we  obtain  the  following  difference  equation 

P(k  + 12)  =  ctP(k  +  11)  +  c2P(k  + 10)  +  ....  +  cuP(k),  (B.4) 
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1) 


P(k)  =  0,  (B.3) 


(Z  - yi^niZ2  - 2cos(ia>Ts)Z 
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where  c.,  with  i  =  are  real  numbers,  which  are  fixed 

given  the  values  of  y,  a,  and  T.  According  to  Eq.  (B.4),  we 
observe  that  P(k  +12)  can  be  expressed  as  a  linear  combination 
of  the  previous  12  measurements  or,  alternatively,  P{k)  can  be 
expressed  as  a  linear  combination  ofP(k-l),  P(k-  2),...,  and  P(k- 
12)  as  follows 


P(k)  =  ctP(k  - 1)  +  c2P(k  -  2)  + ....  +  cnP(k  - 12) 


=  YcP{k-i). 
1=1 


(B.5) 


By  analyzing  the  correspondence  between  Eqs.  (B.5)  and 
(B.2),  we  concluded  that  the  two-process  model  in  Eq.  (5)  is 
equivalent  to  the  12th-order  AR  model  given  by  Eq.  (10),  where 
bi  =  ci,  with  i  =  1,...,12. 
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